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ABSTRACT 

Modern N-body techniques for planetary dynamics are generally based on sym- 
plectic algorithms specially adapted to the Kepler problem. These methods have 
proven very useful in studying planet formation, but typically require the timestep 
for all objects to be set to a small fraction of the orbital period of the innermost body. 
This computational expense can be prohibitive for even moderate particle number for 
many physically interesting scenarios, such as recent models of the formation of hot 
exoplanets, in which the semimajor axis of possible progenitors can vary by orders of 
magnitude. We present new methods which retain most of the benefits of the standard 
symplectic integrators but allow for radial zones with distinct timesteps. These ap- 
proaches should make simulations of planetary accretion with large dynamical range 
tractable. As proof-of-concept we present preliminary science results from an imple- 
mentation of the algorithm as applied to an oligarchic migration scenario for forming 
hot Neptunes. 
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1 INTRODUCTION 

Dynamical timescales in the solar system vary widely. For 
example, Mercury's orbital period is 0.24 yr, Pluto's is 250 
yr, and comets in the Oort cloud can have periods of ~ 10 
Myr, for a difference of over seven orders of magnitude. The 
disparity in timescale increases if we also consider not merely 
orbital period but period at periapse, such as in Sun-grazing 
comets which can come within a few solar radii, or the orbits 
of moons and satellites. The age of our solar system is ~ 
4.6 Gyr, and relevant formation timescales are believed to 
be on the order of a few million years (for the accretion 
of giant planet cores) to tens of millions of years (for the 
formation of the Earth) to hundreds of millions of years (for 
pos sible late-stage rearr angement of the outer solar system, 
e.g. iTsiganis et ai1l2005l h 

The enormous number of orbits required presents a con- 
siderable challenge for numerical studies of planet forma- 
tion. In N-body studies of galaxy formation, by contrast, 
the number of dynamical times is low, and therefore the 
emphasis has been on increasingly complex modelling of 
the physics and ever-larger numbers of particles, both of 
which are amenable to parallelization. Although multipro- 
cessor codes can offer major benefits even for planetary sim- 
ulations, the large number of timesteps limits the particle 
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number to a regime where latency issues loom large. The 
planetary problem is simply, and inescapably, hard. 

The difficulties are yet greater for studies of the forma- 
tion of hot Neptunes, giant planets orbiting very close to 
the parent stars. These planets are unlikely to have formed 
in situ, suggesting that gas disc induced migration will play 
an impor tant role. For exam ple, GJ436b is a hot transiting 
Neptune feutler et al.ll2004h whose density is ~ 1.69 g/cm 3 
jTorredl2007j r This suggests it is a true Neptune analogue 
(i.e. an ice giant, not a small gas giant or a large, rocky super- 
Earth), which raises the question of where its ice originated. 
In the standard models of the minimum-mass solar nebula, 
the snow line beyond which ices can condense is ~ 2.7 AU. 
Therefore one obvious toy scenario for the history of GJ436b 
is that the planet started its life in the outer regions of the 
disc, and then migrated in due to interactions with the gas, 
possibly growing en route. 

At present, the most elegantly constructed general 
integrator for doing lat e-stage N-bod y stud ies of planet 
formation is SyMBA (|Duncan et al.l Il998l h a Kepler- 
adapted symplectic integrator capable of resolving close en- 
co unters, which is de s cende d fro m the original method s 
of IWisdom fc Holmanl i|l99lf ) and iKinoshita et al.l (|l99ll ). 
There are several implementations of the algorithm and its 
variants available, including a parallel versi on which has 
been useful in studying terrestrial accretion (|McNeil et al.l 
2005), and so it would be natural to apply these codes imme- 
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diately. Unfortunately, a direct treatment of the hot Neptune 
problem is completely beyond the reach of standard meth- 
ods, at least at the usual resolution. The orbital period at 
0.05 AU is ~ 0.01 yr, 30 times smaller than at 0.5 AU (a 
respectable inner boundary for studies of the formation of 
the Earth). Since migration over several AU plays a role, we 
cann ot concentrate on a narrow region (e.g. iKokubo fc Idal 
1 19981 ) but must build a more global model, requiring large 
numbers of particles. Furthermore, since the formation and 
migration timescales will be important and often compara- 
ble, some common tricks for speeding up simulations (such 
as increasing the collisional radius) may be dangerous. 

These formation scenarios are of considerable interest, 
and are not easily studied using current numerical and com- 
putational technology; they are messy systems, with many 
non-Hamiltonian forces and events, and not clean celestial 
mechanics problems; and there seems to be a promising di- 
rection for improvement. We therefore seek to develop new 
integrators which will allow us to address these problems. 
Given the difficulties, we are willing to consider approxima- 
tions which we would hesitate to use in other situations, 
such as a detailed study of long-term chaos in the outer so- 
lar system. We will instead sacrifice some precision while 
preserving reliability in the hopes of exploring otherwise in- 
accessible science: in this situation practicality beats purity. 

Fortunately, the very feature which makes the problem 
so challenging - the wide dynamical range - opens up pos- 
sibilities for new methods. SyMBA and its cousins use a 
common timestep for all (non-encountering) objects which 
is set by the minimum pericentric distance. This limitation 
can be worked around when objects only occasio nally enter 
the innermost regions (ILevison fc Duncan 2000) but it re- 
duces to using unacceptably slow Bulirsch-Stoer integration 
when there is always an object in the innermost regions, as 
we expect in our models. However, we note that for an in- 
ner edge of 0.05 AU, a timestep of ~ 0.0005 yr would be 
necessary, but if the inner edge of the problem were 1 AU, 
a timestep of 0.05 yr would suffice. This suggests that us- 
ing the common small timestep results in objects beyond 1 
AU being 'over-integrated' by roughly a factor of 100. If we 
could somehow use the larger timestep for the outer objects, 
then since most particles in simulations of oligarchy tend to 
be in the outer regions (where formation times are longer), 
we might be able to recover something lik e the standard run 
times. Indeed, a decade and a half ago, ISaha fc Tremainel 
(1 19941 ) were already building mixed- variable integrators with 
different timesteps associated with each planet, so there is 
precedent. 

Therefore we set out to construct a new multizone 
method using recent numerical technology which allowed for 
small timesteps in the inner regions and large timesteps in 
the outer regions. We had several desired properties: 

(1) Symplecticity, or at least near-symplecticity, is 
highly desirable. By contrast, time-reversibility (in the ab- 
sence of mergers and dissipational forces) is necessary, as 
many of the good conservation properties of symplectic in- 
tegrators are inherited from their reversibility. 

(2) SyMBA's underlying integrator step is very ro- 
bust, especially at low eccentricities, and has extensive field- 
testing (both in Duncan and Levison's SWIFT and in John 
Chambers's MERCURY). An integrator which reduces to 



this proven approach for objects which are in the same zone 
is preferable. 

(3) Correct close encounter handling is vital. Although 
we may be willing to accept a decrease in encounter accu- 
racy in a few locations (such as at the zone boundaries), 
the vast majority of encounters must be treated using a 
method known to be reliable. More generally, force inac- 
curacies should be limited to distant interactions, as in fast 
force techniques such as treecodes. 

(4) Any discontinuities caused by the existence of dis- 
tinct timestep zones should be kept to a minimum, and at 
or below second order if possible. 

By using existing techniques in the literature, we de- 
velop integrators which meet the above criteria, choosing 
at every branch point the simplest scheme which seems 
likely to work. We combine the Hamiltonian splitting of 
Duncan ct al. ( 1998) with a multistep approach inspired by 
Saha fc Tremainel l| 19941 ). use the transitioning approach of 
Chambers! (|l999l ) to preserve symplectic behaviour, and de- 
rive new transition functions to make the scheme sufficiently 
smooth. 

In iJ5] we briefly review the use of symplectic methods 
in planetary dynamics. In i]2.1l we introduce the Kepler- 
adapted mixed-variable symplectic integrators; in £12.21 and 
M2.3I we explain methods for treating close encounters be- 
tween planets and between a planet and the Sun, respec- 
tively; and in H2.4I we discuss integrators which allow indi- 
vidual timesteps. In H2.5I we construct the new integrators, 
and in £|2.6l we develop appropriate transition functions. We 
presents tests of the method in and a discussion in 
We conclude in §S] 



2 SYMPLECTIC METHODS FOR KEPLERIAN 
POTENTIALS 

Geometric integrators are numerical integration methods 
which attempt to build the properties of the equations 
and their solutions into the integrators themselves, proper- 
ties such a s symmetries a nd their corresponding conserved 
quantities (|Yoshidalll993l ). This matching of geometry be- 
tween the problem and the solver can lead to improvements 
in accuracy and robustness, as well as dramatic increases 
in speed. Symplectic integrators are geometric integrators 
where the geometry of interest is Hamiltonian, and the con- 
served quantity is the natural area element, the Poincare 
2- form dp i A dq 1 , where q and p are the usual generalized 
positions and momenta, respectively, with summation over 
particles i with mass rm.) A symplectic integrator can be 
constructed via operator methods or frequency maps; we 
restrict ourselves to operators. 

For some Hamiltonian H, if we divide it (arbitrarily but 
conveniently) into a kinetic term Ht and a potential term 

Hy, 



H = Ht ~t~ H\ 



(1) 



we can approximate the evolution under H by a second- 
order time-reversible drift-kick-drift 'leapfrog' scheme 



H T ^HZ,' 2 HZ HZ /2 



(2) 



where H£ is the operator generated by the Hamiltonian 
term Hx, and r is the timestep. This integrator solves 
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a nearby 'surrogate' Hamiltonian H — H + H CII where 



0(H C 



t , using 'solves' in the sense of physics and not 



mathematics. (H e rr is a purely formal series which need not 
converge in general, and certainly need not converge at the 
large timesteps used in practice, but is well-approximated 
by its first few terms.) 



2.1 Wisdom-Holman method 



The w ork of IWisdom fc Holmanl (|l99ll ) and lKinoshita et all 
l|l99ll ) sparked a dramatic revolution in planetary integra- 
tions by specializing the general techniques of symplectic 
integration to the unique properties of Keplerian dynamics 
in the solar system. It is well known that the orbits of a 
system of small bodies around a large central mass are very 
nearly conic sections, and this Keplerian motion can be ad- 
vanced (almost) analytically. This suggests that instead of 
dividing the Hamiltonian into kinetic and potential terms 
as in the standard leapfrog, we should divide it into a Kep- 
lerian term and a term corresponding to the perturbations 
between the planets: 



H — Hkc P + Hiv 



(3) 



The difficulty arises in finding a canonical set of coordi- 
nates in which the N-bod y Hamiltonian takes this shape. 
IWisdom fe Holmanl (|l99ll ) discovered that using Jacobi co- 
ordinates succeeds. In this coordinate system, the position 
and momentum of the jth object (of N planets, where the 
masses are given by the rrij , and the Sun is j = 0) are given 
relative to the centre of mass of the inner bodies, i.e. those 
with index k < j. They write 



H. 



Kep 



E 



G m-j mo 



2m' 



and 



E 



/ Gmjino Grrijmo 



9jo 



(4) 



JV-l JV 

ST Gm i rnk (5) 



3=1 k=j + l 



9jk 



where the primed quantities are Jacobi coordinates and 
lik — 1 9 j — IkY in tne absence of close encounters, _ffi nt <C 
Huep and the non-Keplerian perturbations on the orbits 
due to mutual interactions are small. We can construct 
a second-order integrator as before, where the new Hkc P 
drifts are 'rolls' along the Keplerian conic section, and 
the Hi n t kicks are the Cartesian perturbations between the 
planets; hence the nam e 'mixed- variable symplectic' (MVS; 
ISaha fe Tremainei ri992) has been used for integrators of this 
type, as the integrations are in effect carried through in 
both Cartesian and Keplerian variables. This approxima- 
tion to the orbit is vastly superior to the linear-path ki- 
netic/potential decomposition, and the symplecticity pro- 
vides robustness and stability. As a result, one can use much 
larger timesteps than would be otherwise permissible, mak- 
ing planetary simul ations running for t he age of the solar 
system feasible (e.g. iDuncan et al.il 19951) . 

For a brief review of the history of various map- 
ping methods in s olar s ystem dynamic s, see ch. 9 of 
iMurrav fe Dermottl (|l999r i. IWisdoml (|2006t ) also provides a 
useful review. 



2.2 Close encounters between planets 

Despite its many benefits for studies of well-separated plan- 
ets for long timescales, the use of Jacobi coordinates in the 
MVS approach limits its applicability. Studies of planetes- 
imal accretion, a highly chaotic and stochastic process, re- 
quire the ability to handle radial reordering of objects and to 
resolve close encounters. In the Wisdom-Holman mapping, 
the coordinate frame depends upon a fixed ordering of the 
objects, and cannot be updated as the system changes with- 
out breaking symplecticity. There is also no natural way to 
treat encounters: first, by varying the timestep you lose the 
symplecticity of the integration as the composition of two 
symplectic steps need not be symplectic; and second, even 
if one could vary the timestep, the coupling between objects 
in the Jacobi scheme means that all external objects are 
also implicitly involved, so an encounter cannot be treated 
independently. We would prefer an integrator in which each 
planet was treated equivalently and which could adapt the 

timestep when planets strongly interact. 

After some early experiments (e.g. iLevison fc Duncan] 
1994) . a super i or sol ution to the problem was presented in 
Duncan et al.l (|l998h (hereafter DLL98) which used an in- 
genious choice of coordinates to avoid the coupling between 
planets, a particular choice of Hamiltonian splitting to en- 
sure that all planets were treated equally, and a clever de- 
composition of the potential into shells to allow effective 
changes in the timestep. The authors developed a 'demo- 
cratic heliocentric' (DH) method , where 'democratic' means 
'symmetric with respect to the labelling of the planets', and 
which avoids the particle entanglement of Jacobi coordi- 
nates. In these coordinates, one uses heliocentric positions 
but barycentric momen ta. (Note that t he same system is 
called 'mixed- centre' b y [Chambers 1999 and 'canonical he- 
liocentric' by IWisdoml ETooe.fl To be explicit, from q and 
p they compute new conjugate coordinates Q and P, such 
that 



i 



EN 



i = 
i / 



(6) 



and 



This results in a new form for the Hamiltonian, 

H{Q 1 ,P l )= tfsun + #Kcp + H int (8) 

where 



1 

2mo 



E^ 



\Pi\ Gm,m 



Irrii 



! = i 

n — 1 n 



\Qi 



x ^ x ^ Grrnrxij 



(9) 
(10) 

(11) 



1 The temptation to suggest yet another name is resisted only 
with difficulty. 
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-Hsun generates a linear drift of particles' positions (tak- 
ing its name from the fact it is determined by the barycentric 
momentum of the Sun) , Hkc P corresponds to a pure Kepler 
orbit, and -Hint, as before, is the term due to the interaction 
of the particles. Note that - unlike with the Wisdom-Holman 
mapping (eqs. [4] and [5]) - an encounter between two parti- 
cles is separable, in that the terms involving the positions 
and momenta of the two objects can be pulled out of Hk bp 
and Hi n t. (In contrast, -ffs un is not separable in this way; we 
return to this subject.) 

As IWisdoml (2006) notes, the same canonical he- 
liocentric coordinate system was used in previous work 
ijTouma fc Wisdom|[l993l.ll994h. with differ e nt Ha miltonian 
splittings. For example. iTouma fe Wisdom! (|l993l ) used the 
alternate splitting 



1 n 

#Sun = — V PiPj 



i<i<j 



2 m 



Grriimo 



(12) 



(13) 



where p is the reduced mass. This splitting has the advan- 
tage that it preserves Kepler's semimajor axis-period rela- 
tionship. However, the use of the reduced mass means that 
two objects with the same position and velocity but different 
masses will experience different Kepler drifts and different 
linear drifts, which can present difficulties for treating close 
encounters (DLL98, §4). In this sense the DLL98 splitting is 
more 'democratic', at the cost of being slightly less accurate 
for a given orbit. Accordingly, we will use 'canonical helio- 
centric' for the coordinate system itself, and reserve 'demo- 
cratic heliocentric' for the specific three-term Hamiltonian 
splitting used in DLL98 and described by eq. [5] 

From this DH splitting, they construct a second-order 
integrator which will be the basic DH step. Following the 
conventional notation, let L be the operator generated by 
Hsun (L for 'linear drift'), let K correspond to H lrA (K for 
'kick'), and let D correspond to Hkc P (D for 'drift'). Using 
this form, the system is advanced one timestep r by applying 



H T ^L T ' 2 K T/2 D T K T/2 L T/2 



(14) 



As naive adaptive timestepping changes the surro- 
gate Hamiltonian and thus breaks symplecticity, DLL98 
develop a technique involving recursive subdivision of the 
timcsteps along with decomposing the force into a se- 
ries of shells around the particles, and associating each 
force shell with a different timestep. Similar approaches 
had been attempted pr eviously in molecular dynamics 
i|Skeel fc Biesiadecki[i994 ) with limited success, but DLL98 
realized that there were relevant smoothness constraints on 
the transition from shell to shell. 

Starting from the basic DH step, we relabel K, D as 
Kq, Dq. We replace the D T by a sequence of M, smaller sub- 
steps, \K T J 2Ui D T J Ui Kj^ 2 4 ] M< where Ki is an interaction 
term to be defined, and then repeat the process indefinitely 



at higher index: 

H « U' 2 Kl 12 Dl Kl 12 U' 2 

« U' 2 K T /2 [Kl /2Ml Dl /Ml Kl /2Ml ] M i Kl /2 U' 2 

~ r/ 2 k t/2 \k t/2Mi 

l^r/ajfilfj jjt/MxMz ^ T /2MiAf 2 jM 2 > 

K T/iM 1]Ml K T/l LT/ 2 



Each of the above integrators, as well as the i — > oo limit, 
has a fixed surrogate Hamiltonian; the timestep r never 
changes. Since the Di terms commute with themselves, then 
if succeeding Ki terms are merely the identity operators the 
neighbouring Di terms reduce to one large drift. For exam- 
ple, if Ki = for all i > 0, then the integrator reduces 
to the basic DH step with timestep r. However, if K2 were 
nonzero, then since Di and Ki do not commute, the re- 
duction would not occur and the integrator would act on 
the smaller timestep. In this scheme, an infinite number of 
terms are always active, but need not actually be considered 
unless the intervening Ki terms are nonzero. The problem 
becomes finding a decomposition of the potential H lnt gen- 
crating the K terms Kq, Ki , K2, ■ ■ ■ which will (1) recover 
the standard DH step when no close encounters are occur- 
ring and the larger timestep r suffices, i.e. Ko — K, and 
Ki — for i > 0; (2) always sum to the correct amount of 
force, i.e. Ki = K; and (3) do so in a sufficiently smooth 
fashion. To construct such a decomposition, DLL98 imagine 
a decreasing sequence of shell radii Ri > R2 > . . . with a 
corresponding decomposition of the potential into V\, V2, • • ■ 
where the magnitude of the potential terms Vk smoothly 
decreases to zero outside the shell range. By studying the 
error Hamiltonian they find several necessary conditions on 
the decomposition relating to the smoothness of the transi- 
tion of the potential from term to term. 

The resulting algorithm SyMBA, combining the DH co- 
ordinate system with the recursively subdivided smoothly- 
transitioning zone operators, works remarkably well. Never- 
theless, it is challenging to implement and difficult to test. 
For example, a correct version of the algorithm will often 
show worse energy and angular momentum conservation on 
a given encounter than a subtly incorrect one, requiring large 
test suites; and since the solution is very finely tuned to the 
pla netary prob le m, of f-the-shelf routines are of little use. 

IChambersl (Il999l ) developed a much simpler though 
more expensive approach. Of the three major advances in 
SyMBA - DH coordinates, the smoothness conditions on 
motion between operators, and the recursive subdivision of 
the Hamiltonian - only the first two are strictly necessary 
to the resulting algorithm. The major benefit of the recur- 
sion is that it pushes the transitions entirely into the kick 
operators and preserves the ability to apply the drift oper- 
ator merely by solving Kepler's equation, albeit at the cost 
of considerable complexity. If we surrender this requirement, 
however, then we can decompose the Hamiltonian into three 
terms 

H = H Sun + [F(i,j) Hij] + [(1 - F(i,j)) Htj + H Kcp ] (16) 

where Hij is the component of H ln t involving objects i and 
j, with implied summation over all i,j with i < j. F(i,j) 
is a transition function which is 1 when the objects are dis- 
tant and approaches when they are undergoing an en- 
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counter. As in the basic DH step (eq. 114) 1. this decompo- 
sition gives rise to a five-operator step, but these are the 
only five terms to consider, unlike the much larger number 
of potentially active terms in eq. 1151 It is true that when- 
ever an encounter is occurring and F(i,j) could be nonzero 
then the interaction terms can no longer be advanced ana- 
lytically, but they can still be numerically integrated to high 
precision using standard techniques such as Bulirsch-Stoer. 
Moreover, and this advantage should not be underestimated, 
testing that the above algorithm is implemented correctly is 
far more straightforward than testing SyMBA. Whenever 
the cost of the numerical integration due to encounters is a 
small fraction of the computation, the Chambers-style ap- 
proach (although perhaps not as elegant as SyMBA itself) 
is likely preferable on pragmatic grounds. It should also be 
noted that despite the use of numerical integration to han- 
dle the transition-weighted terms in the Hamiltonian, the 
mapping itself remains symplectic; or as symplectic as any 
float ing-point im plementation of sufficient precision can be. 
rSee lSkeelll999l for a general technique to recover symplectic- 
ity when using otherwise non-symplectic approximations.) 



2.3 Close encounters with the Sun 

In any case, both SyMBA and the Chambers variant have 
been successfully applied to many studies of the later stages 
of planet formation. Their chief weakness is that neither 
can easily deal with objects undergoing close encounters not 
with each other but with the Sun, such as high-eccentricity 
Sun-grazing comets. Unlike the case of mutual encounters 
in which the important terms are -Hint and Hkc P , during a 
close solar passage ffsun must be evaluated more frequently. 
Therefore, to resolve such orbits correctly, one must choose a 
timestep small enough to resolve the pericentre passage, and 
that timestep must be fixed for the entire integration, even if 
such encounters are very rare. To overcome this limitation, 
iLevison fc Duncan! (|2000l ) added a Chambers-style splitting 
on top of SyMBA (here we suppress the planetary encounter 
terms, which are handled as described before) and use 

ir = (i-f)h;£+ 

(F Hsun + H Kep ) T + (17) 

H = (1 - F) H S un + Hint + (F H S un + #Kcp) (18) 

where F is a transition function which is when no object 
is near the Sun and 1 when any object is. In this scheme, 
when any object is undergoing a close solar approach, then 
the work of performing an integration step is pushed into 
the new Kepler step which is must be handled numerically. 
The non-separability of -ffsun - i.e. the inability to isolate an 
individual object as in i/xep - becomes a serious inconve- 
nience here, as it requires the numerical integration of every 
object in the system even if only one o bject enters the in- 
ner zo ne. Nevertheless, as explained by ILevison fc Duncan! 
(2000), this approach provides the speed of SyMBA when- 
ever the inner regions are empty and yet can successfully sur- 
vive occasional interlopers (such as objects due to be ejected 
during violent periods in the formation process). 



This method cannot be directly applied as a solution for 
the numerical challenge of hot exoplanet formation. If we set 
the inner boundary at a typical terrestrial-formation value 
like 0.5 AU, then there will usually be many, and almost al- 
ways be some, protoplanets and planetesimals in the inner- 
most region. This means that every particle's Hkc P will be 
numerically integrated on every step, and the speed benefits 
are lost. Moreover, if there are mutually-gravitating objects 
inside the inner boundary, then the situation is worse: the 
above integrator will evaluate the interactions between bod- 
ies not undergoing close encounters on the outer timestep r, 
which may bear little relation to the dynamical times for the 
inner objects. This particular problem could be corrected by 
bringing Hint under the transition function F as well, but 
then we have merely recovered - in an impressively round- 
about fashion - a Bulirsch-Stoer integrator. 

2.4 Individual timesteps 

Recognizing that in an MVS integration of the solar system, 
one is taking hundreds of times more steps per Pluto orbit 
than would be necessary if not for th e presence of the inte- 
rior planets, ISaha fc Tremainel l| 19941 ) construct a leapfrog 
integrator with individual timesteps. They split Hkc P and 
-Hint into N terms each, such that 

JV JV 
Hkc p = HKc P ,i , H ln t = ^^#int,i (19) 

i = l i = l 

where HKc P ,i is the Kepler term for planet i (with objects 
labelled in increasing semimajor axis), and -ffint.i is the in- 
teraction term between planet i and planets i + 1 through 
N. Translating into our notation, we use Di and Ki to refer 
to the operators as before, and assign a timestep t; to each 
planet, where the largest timestep t = tn and Ti+i/Ti = Mi 
for Mi an integer. Starting with the Hamiltonian only in- 
volving the outermost planet (note that Kn is the identity 
operator): 

H T » D t n n/2 K t n n D t n n/2 (20) 
and recursively applying 

KI -> [D^K^KID^Y- 1 (21) 

for i > 1, one obtains an N-level integrator. More concretely, 
consider a two-planet case with timestep ratios of 1:3. The 
resulting integrator (after removing the Ki term which does 
nothing) 

H T « D T 2 2/2 [Dl l/2 KI 1 Dl l/2 f D T 2 2/2 (22) 

has a timestep of n for the inner planet, and T2 for the outer 
planet, as desired. Note that the integrator is time-reversible 
even though the objects are not synchronized with respect 
to D when the interaction term K\ is applied. (Also note 
that in practice one would combine neighbouring D\ terms.) 

It is important to recognize that the timesteps associ- 
ated in this method are individual but not adaptive; they 
must be set at the start of the integration, and attach not 
to spatial zones but directly to objects. It is therefore un- 
able to handle migrating objects, and as an MVS method 
inherits the previously mentioned weaknesses of Jacobi co- 
ordinates for our purposes. However, it demonstrates that 
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asynchronous multi-stage integration in MVS-like contexts 
can be constructed. 



2.5 Constructing the new integrator 

We now have the necessary ingredients to construct 
a symplectic integrator, Naoko ("New Adaptive Or- 
thochronous Kepler Orbiter"), which is Kepler- adapted, 
close-encountering, and yet allows for zones with different 
timesteps. Recall the basic DH step: 



H T ^L T/2 K T/2 D T K T/2 L T/2 



(23) 



We will seek a generalization of this step for the 
multiple-zone case. We will define our timestep zones by 
dividing the system into radial shells such that an object's 
instantaneous heliocentric radius determines its zone assign- 
ment. (This definition is mentioned here for concreteness; 
other choices are possible. We motivate this particular choice 
in £12.61 ) We label the zones using integer indices, starting 
at for the innermost zone, and an integer subscript on 
an operator corresponds to the operator for that zone, in 
a sense which will be made explicit later. For example, Dq 
advances all objects in zone for a timestep r under Hkc P 
but does nothing to objects in other zones: the commuta- 
tor bracket \Dj,Dj] = 0. Note that this differs from the 
ISaha fc Tremainei i|l994h usage of subscripts to refer to plan- 
ets. 

To begin with, we defer consideration of interactions 
between the planets. Under this simplification, and starting 
with zone 1, eg. 1231 becomes 



H T « L T/2 Dl L 



-t/2 



(24) 



The L terms cannot be divided into zones as Hsun is not sep- 
arable. Nevertheless, we can incorporate more zones by ap- 
plying the individual-leapfrog approach of lSaha fc Tremainei 
(j 19941) and recursively subdividing the L steps using this ex- 
pression. That is, we can write 



L r/2 

resulting in 



L T/i Dl /2 L T/A 



M 



^T/4Jf jjt/2M £ T /4M j" 

^t/«„ Dl /2M ° L T/4M °j M ° 



(25) 



(26) 



where Mo is an integer setting the number of zone steps 
per zone 1 step. 

This integrator looks promising. In the absence of any 
objects in zone 1, this is merely 2Mo (kick-free) DH steps 
of size t/2Mo next to each other, and in the absence of any 
objects in zone 0, the Do operators do nothing and the L 
operators collapse, reducing to a DH step of size r. There are 
2Mo zone drifts per zone 1 drift, so if we assume each drift 
takes equal time then if Ni > 2MqNq, the zone 1 compu- 
tation dominates. As long as Ni 3> 2MqNq, then handling 
objects in the innermost zone - far from requiring a ma- 
jor decrease in system timestep as in the standard approach 
- is effectively free; and, importantly, they can be handled 
independently of the outer objects. 

Despite appearances, the presence of the L terms on 
the smallest timescale does not remove this separability. 



Although L must be formally applied with the innermost 
timestep, this does not mean that L must be computed at 
that frequency. Since Do does not affect objects in zone 1, 
and L changes their positions but not their velocities, we 
need only determine their contribution to L at the begin- 
ning of the zone substep, and we can delay actually moving 
them until the beginning of the zone 1 substep. Implement- 
ing a lazy-evaluation scheme is relatively simple. (Since L is 
so trivial to evaluate and apply, it is seldom a bottleneck, at 
least in the serial case. For parallel implementations, we find 
that lazy evaluation is vital, because otherwise the commu- 
nication overhead involved produces enormous scaling diffi- 
culties. Potential implementors should bear this warning in 
mind when designing data structures.) 

One can easily generalize to more zones. The three-zone 
case is the simplest integrator where the transition functions 
we will develop can be carried through to a larger number of 
zones, and therefore we will use it as our standard example. 
In a minor abuse of notation, define S T = L T ^ 2 D T L T ^ 2 (note 
that unlike with L, D, and K, [S T ] M S Mt .) 



Then we have the two-zone integrator 



and the three-zone version 



H T = 



/4M M! 



M n Mi 



Dl 



[[«' 



T/±M Mt] M ° r ,T/2M 1 \ g t/4,MoMi-\ M °] Mi 



(27) 



(28) 



[So* 



where Mi sets the number of zone i substeps per zone i + 1 
step. 

After one step, objects in all zones have been advanced 
the correct total time under each operator. Moreover, ne- 
glecting the influence of objects in other zones, each zone 
has experienced what is locally a standard DH step: objects 
in zone took a timestep of t/4MoMi; objects in zone 1 had 
a timestep of t/2M\\ and objects in zone 2 had a timestep 
of T. 

We must now choose where to place the force operators. 
To simplify the discussion, we set Mo = Mi = 1, and assume 
that the zones are separated such that r, t/2, and r/4 are 
appropriate DH timesteps for all objects in the respective 
zones; recovering the general case is straightforward. The 
kick-free three-zone integrator is given by 

H T = I/l*Dl IA Ul*D\ l2 Ul*Dl IA I/l*Dl 

u/&D r /iu/&D r /2u/&D r /iu/& (29) 

Let Kij be the kick operator between zones i and j. In 
order to treat close encounters using the Chambers splitting, 
we must have each Di surrounded by two Ku. 

H r = L^ S K^ 8 D T /4 K^ S L^ S 
Kl( 4 Dl /2 Kl( 4 

r-r/8 K r/8 d t/4 k t/8 , t/8 
1j ^00 ^o "oo lj 

Kli 2 DlKU 2 (30) 

This integrator does not incorporate interzone forces, but 
intrazone forces are evaluated on the appropriate timescale 
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relative to the drift timescale for objects in that zone. One 
can therefore apply the Chambers-style encounter handling 
between the sets of KaDiKa as described in section 12,21 
and we recover the standard approach. (In some situations 
we have found the DKD operator splitting is superior, but 
for reasons involving implementation details we will restrict 
ourselves to discussing the KDK version.) 

There are many possible arrangements for force commu- 
nication between zones. We optimistically choose the scheme 
requiring the fewest force calculations, and update forces be- 
tween zones i and j on the outermost zone's timestep. This 
choice results in 

H T ~ L t ' 8 K^ 8 D t /4 kXl t/8 

„t/4 „t/4 n r/2 „t/4 „t/4 
T t/S r t/S D T/i K T/8 j-r/8 
/ 2 If*/ 2 r>r IfT/2 K T/2 



-"02 A i: 



2 K22 DIK^ 2 



K, 



72 



(31) 



2-^22 xl 12 
tt/8 k t/8 d t/4 k t/%jt/s 

u -"-00 ^0 -"00 u 

k t/4 n r/2 k t/4 jst/4 

"01 A xi u \ -"11 -"-oi 

TT/& js-t/8 J-.T/4 k t/8 tt/8 
L A 00 U -"00 ^ 

This approach succeeds in separating the integration of 
zone i objects from zone j objects. Objects in zone 2 need 
only have drifts and kicks involving them evaluated on the 
timestep r; objects in zone 1 on r/2; and zone on r/4; 
and as already discussed L is not a problem. Note that the 
above integrator was built to minimize the number of force 
operators, but one could evaluate the cross-zone kicks [K\2, 
for example) more frequently if desired. 

The above integration technique obeys Newton's third 
law regarding interparticle forces. Although forces between 
objects in different zones are computed less frequently than 
forces within a zone, whenever gravitational accelerations 
are computed between two bodies the accelerations are equal 
and opposite (even for close-encountering objects in a tran- 
sition zone, to be discussed later) and they are immediately 
applied and turned into changes in velocity. No force lag or 
accumulation is involved. 



2.6 Transition functions 

We recall that the key insight of DLL98 is that to preserve 
symplectic behaviour while effectively changing the integra- 
tion step (and therefore the surrogate Hamiltonian) an ob- 
ject must experience a smooth transition from one integra- 
tion regime to another. Roughly speaking, if an object's 
transition is sufficiently smooth, then instead of suddenly 
finding itself evolving under a different surrogate Hamilto- 
nian, it believes it is merely in a different regime of the 
original Hamiltonian. Accordingly, we now return to the 
previously-deferred issue of choosing appropriate transition 
functions, which turns out to be the most challenging part 
of the problem. 

Let f(x) be a switch function, a real-valued nondecreas- 
ing function defined on [0, 1] with /(0) = and /(l) = 1. 
We take the extension outside this domain (equal to below 
and 1 above) as given, under which convention f(x) — x is 
a switch. The simplest switch is a shifted step function: 



/step (as) = 



x < 1 
x> 1 



(32) 




Figure 1. Various candidate switch functions, discussed in 32,61 

Hamiltonian from one term to another, and symplecticity is 
lost. DLL98 suggest 

/dll 3 (z) = 3x 2 - 2x 3 (33) 

as one of their switches, which has f'(x) — at both end- 
points, as well as the higher-order 

3\ 



/dlltOe) = £ 4 (35 - 84a; + 70x 
IChambersl (Il999h suggests 
fch{x) =x 2 /(2x 2 -2x + l) 



20ar 



(34) 



(35) 



as a useful co mpromise between smoo thness and efficiency 
of evaluation. iRauch fc Hoim an (1999) prefer 



/rh(x) = - ( 1 +tanh 



2x 



(36) 



x(l — x) 

for which all derivatives vanish at the endpoints. Figure [T] 
shows the various functions. We will use the cubic polyno- 
mial switch /dll3, but the construction is independent of 
this choice. (This issue is discussed further in M3.1.1H 

To reduce clutter we define a rescaling function on the 
switch, 



C(x, xo, xi) = /dll 



x 



Xi - x 



(37) 



However, this leads to sudden movement of portions of the 



and build our transition functions from this base. 

As described briefly in £|2 . 5 1 we divide the system into 
zones based on instantaneous heliocentric radius r, and 
imagine a set of spherical shells around the Sun. Each zone 
has a transition region near the interior and exterior edge of 
the shell in which operators involving both adjacent zones 
will act on an object. Let FU be the locations of the i zone 
boundaries and 2hi be the widths of the transition regions, 
as illustrated in figure [2] 

An unfortunate consequence of this choice is that the 
resulting evolution of transiting objects under D is no longer 
analytically integrable. It is tempting to construct a tran- 
sition function based not on r but on the osculating semi- 
major axis a, which is constant during D and can therefore 
be used (along with all other orbital elements save those 
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W D (r, 0) 
W D (r, 1) 
W D (r, 2) 



Ro 



i i r~ 

R,-h, R, R, + h 



I I I 
■i R2 — h 2 R2 ^2 ^2 

heliocentric distance 



Wn(r,i) = < 





C(r,Rr,Rl 
1 

1 - C(r, 




R~ ^r^ R 
Rt -R 



t+i 



(38) 



Rt+i < r 



Figure 2. A sketch of the radial zone scheme. 



where r is the instantaneous heliocentric radius and i is the 
zone index. This is nothing more than a SyMBA-style tran- 
sition applied to r instead of to the interplanet separation. 
Although in general each Ri can have an associated transi- 
tion zone, in practice the boundary 'transitions' are slightly 
degenerate. Typically one would set the inner transition re- 
gion, (Ro, ho), and the outer transition region, here (R3, /13), 
well outside the radii of interest so that the sum of weights 
over all zones is equal to 1 for all integrated objects regard- 
less of r. In practice we make Ro smaller than our inner 
edge (set by the physics of the problem or the numerics of 
our timestep), R3 larger than the outer edge, and let ho and 
/13 be some arbitrary small distance. Figure [2] sketches the 
resulting scheme, with three timestep zones and effectively 
two (not four) transition zones. 



such as the mean anomaly describing the position along the 
orbit) to build a weighted but integrable D. Indeed, exper- 
iments show that an a-based function can work for isolated 
objects in the transition region. However, doing so vastly 
complicates the treatment of close encounters, as two en- 
countering objects must have similar instantaneous r but 
can have very different a: consider a high-e zone object 
at apocentre meeting a high-e zone 1 object at pericentre. 
Under an r-based weighting scheme, the difference in their 
effective drift timestep is bounded by their physical radial 
separation which gets smaller as the encounter gets deeper, 
whereas under an a-based scheme the drift timestep experi- 
enced by the inner object could remain very different from 
that of the outer object, which is not a recipe for numeri- 
cal stability. A function which smoothly changes dependence 
from r to a with changes in the encounter status may be pos- 
sible, but in our view the potential benefits are outweighed 
by the resulting complexity. 

For the integrators presented here, such as that of 
eq. 1311 there are three types of transitions we must con- 
sider: (1) those involving drifts; (2) those involving distant 
kicks; and (3) those involving close encounters. Each has its 
own peculiarities, and treating (2) and (3) simultaneously is 
rather awkward. 



2.6.1 Drift transitions 

The drift transitions are the easiest to handle. We need a 
function which will yield (for example) the full D\ term for 
an object completely within zone 1, and likewise with D2 
and an object in zone 2, but which will return smoothly- 
varying intermediate values for objects inside the transition 
zones [R\ — hi, Ri + hi] and [R2 — hi,R2 + ha], and gen- 
erally [R~ , Rf] where Rf 1 = Ri ± hi. (For simplicity we 
will restrict ourselves to considering symmetric functions, 
although since objects in inner zones will have higher ve- 
locities and different effective timesteps, it is possible that 
an asymmetric function could yield better results.) Such a 
function is given by 



2. 6. 2 Distant forces 

Treating the forces is more difficult, as it involves not only 
the interacting bodies' two distinct orbital radii but the sep- 
arations between them. We seek a function which will ensure 
the objects experience sufficient continuity in their forces as 
they move from zone to zone. We will first consider only 
distant forces (i.e. we imagine all forces are 'soft', requir- 
ing no special attention) and then correct to handle close 
encounters. 

Consider an integrator with three zones, labelled i, j, 
and k from innermost to outermost, and two planets, with 
orbital radii ro and n. Let the objects start in the zone j. 
Initially, all the weight should be in Kjj . If we keep ro = ri 
and move the pair together both inwards and outwards, we 
see that we need smooth transitions between Ku and Kjj as 
well as between Kjj and Kkk- If we instead keep ro constant 
while increasing n, then as the second object crosses from 
j to k, then the Kjj weight must decrease and Kjk must 
increase until ri is fully within k and the Kjk weight is 1. If 
we decrease ro and increase ri, then when both objects are 
in the centre of the transition zones (i-j and j-k), we will 
need to spread the force out over four different terms: Kij, 
Kj k , K ih , and Kjj. 

This can be achieved by symmetrizing Wd' 



W diBt {ro,ri,i,j) 



(W D (r ,i)Wn(ri,j) + 
W D (n,i)W D (ro,j))/2 

W D (r ,i)W D (ri,j) + 
W D (ri,i) W D (ro,j) 



(39) 



2.6.3 Close encounters 

Now we must ensure that close encounters are correctly 
shared. This requires a transition function W c \ which will 
move the force between two objects from the Ku term to 
the Di term as in the Chambers methods, resulting in a 
triplet of terms 



\w A K v 



+ [(1 - W c i)Ku + Di 



\w cl Ku 



(40) 
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(eq. I45p to the multilevel step of eq. 1311 The order of the 
substeps and the relationships between the r; are deter- 
mined by the recursive subdivision of eq. [5T] (as in eq. 131 [) . 
For each substep of level I, the algorithm proceeds as follows: 

(1) Do the linear drift: if I = 0, apply Lo for To/2. 

(2) Apply the distant forces between I and all interior levels: 
apply Wx(«, I) Wd Ku for all i ^ I for n/2. 

(3) Drift (including possible encounters): apply the 
distant forces between I and all interior levels: apply 
Wk(1,1)(1 - Wcl)Ku + Di for n. (Note that any imple- 
mentation would avoid doing this numerically when the 
kick term was known to be zero because the objects are too 
well-separated.) 

(4) Apply the distant forces between I and all interior levels: 
apply I) Wd Ku for all i ^ I for n/2. 

(5) Do the linear drift: if I = 0, apply Lo for ro/2. 

A moment's consideration confirms that as promised 
in j{T] this is not h ing m ore than the step subdivision of 
ISaha fc Tremainel |l994) applied to the Hami ltonian of 
DLL9 8 with the close encounter handling of I Chambers] 
(1999), with staggered force calculations, and some new 
transition functions incorporated to ensure smoothness. 

In practice, we do not use the above scheme directly, 
but instead use a 'lower -level' scheme (motivated in part by 
iRauch fc Holmani ri999') in which we do not apply the tran- 
sitions between terms in the Hamiltonian but between the 
resulting dq/dt and dp/dt. That is, given two pieces of the 
Hamiltonian Ha and Hb and a transition function / = f(q) 
between them, instead of starting with the decomposition 

[fH A ] + [(I- f)H A + H B ] (46) 



eq. \W\ We take 

W cl = C(ds,ds clit /2, ds crit ) (41) 

where ds cr it, after DLL98, is several times the sum of the 
Hill radii (with a possible additional dependence on orbital 
velocity; as Chambers has noted, what matters is ensuring 
you integrate through the transition .) Note tha t unli ke the 
transition function recommended in IChambersI dl999l ). W cl 
moves all of the encounter into the numerically-integrated 
term (1 — W c \)Ku + Di at separation ds — ds cr it/2, not at 
ds = 0. 

2.6.4 Combining the transitions 

Each of the three above transition functions - Wd, which 
applies to the D operators; Wdist, which applies to all Kij 
operators; and W c i, which applies between Ku operators 
- makes sense independently. The natural way to combine 
them would give (for example, suppressing indices on the 
transition functions): 

(W/ dist Af23) T3/2 (W dist W cl K 33 ) T ^ 2 

[W dist (l-W cl )K 33 + W D D 33 ] T3 (42) 
(W dist W cl K i3 ) T ^ 2 (W diBt K 23 ) T ^ 2 

The above approach has a minor problem, however. 
Consider two objects undergoing an encounter in a tran- 
sition zone (say, the i-j boundary.) The above functions will 
attempt to share force across three terms: Ku, Kjj, and the 
cross-term Kij . However, the integrator is only built to treat 
the encounter using the KuDiKa and KjjDjKjj substeps. 
Any force that the Kij operator is assigned is only sampled 
on the larger timestep, and if the two objects are both in 
the middle of the transition zone this can be as much as 
half the total force. This will result in a highly inaccurate 
integration. 

How can this be repaired? In the case of an encounter, 
Ku and Kjj should share all the force: 

f (W D (r ,i)W D (r 1 ,j)+ i = j 
WW(r ,ri,t,j) = ^ W D {r u i)W D (r ,j))/2 (43) 
{ i + j 

In the absence of encounters, the above method (eq.!42p 
should work. Therefore we define yet another transition 
function, 

W4hift = C(ds, dsc-it, 2 ds crit ) (44) 

and use it to smoothly interpolate between the no-encounter 
case when all operators are involved and the encounter case 
when only the intrazone operators share the force. (Here 
we will require that objects cannot undergo mutual encoun- 
ters unless they are in neighbouring zones. Removing this 
limitation is possible but unnecessary for our intended ap- 
plications.) Combining the above, we write 

W K (i,j) = WWftWdist + (1 - HWfOWVar (45) 

and replace instances of Wdist with this encounter-corrected 
expression. 

2.7 Assembling the integrator 

The final integrator is constructed by applying the drift 
transition function (eq. 138 p and the kick transition function 



which produces 

d l = UfH A ) , f| = -|-(/tf A ) 

d l = ~ f)HA +H B ) , f = -&((1 - f)H A + H B ) 

as the derivatives to be integrated, we write 

dq r dHj\ dp r 9^4 

di ~ dp ' dt ~ J dq (A 

^ = (l~f)^+H B ,f t =-{l-f) 9 -§f+H B { 

and apply the transition directly between the derivatives 
themselves. This saves some computation (of the cross-terms 
Hdf/dq, at least; the Hdf/dp terms are all because / 
is not a function of the momenta) and increases the effec- 
tive smoothness of the transition for a fixed / because one 
does not lose a degree of smoothness passing from H to its 
derivatives. It might be objected that this means the result- 
ing integration technique is no longer strictly symplectic, 
but in fact one can construct a Hamiltonian (admittedly 
somewhat artificial) for which an obviously symplectic in- 
tegrator generates exactly this algorithm. However, even if 
this were not the case and we were to view the integra- 
tor as merely 'near-symplectic', it remains time-reversible 
and works well in practi ce, which agrees with the results of 
IRauch fc Holmanl |l999) who applied the transition func- 
tions in their SyMBA-style decomposition to the forces and 
not the potential. Of course, one may work with the original 
Hamiltonian-level splitting if preferred. 

One easily overlooked issue which is only apparent when 
considering a global simulation is that in the SyMBA scheme 
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the (effective) Hill radii - the Hill radii used in the close en- 
counter criteria - are fixed at the start of the integration. 
This is necessary to ensure that the integration remains both 
symplectic and analytically soluble. However, in cases where 
there is considerable inward migration, an object will have 
an unnecessarily large close encounter criterion when it ar- 
rives at the inner portion of the disc. Similarly, an outward 
migrating object could miss encounters. The strict solution 
is to bring the Hill radius dependence of the close encounter 
weight function under the numerical integration; the lazy 
not-quite-symplectic solution is to update the effective Hill 
radii on some interval, at the risk of interfering with ongo- 
ing encounters. All things being equal, if the system is messy 
and the number of updates is low it is unlikely to make a 
statistically significant difference. 



3 TESTS OF THE METHOD 

The algorithm w as implement ed using the existing code- 
base of miranda (|McNeill l2006) as a framework, and using 
a standard Bulirsch-Stoer routine to handle the numerically 
integrated terms in the Hamiltonian. Since the method re- 
duces to the well-understood Chambers approach in the one- 
zone limit, we will concentrate on testing the performance 
of the multiple-zone aspects of the code. We have confirmed 
by comparison with miranda in both SyMBA and Bulirsch- 
Stoer modes that the new code behaves as expected in the 
one-zone case. Unless otherwise specified, the numerical in- 
tegration tolerances were set at 10 -16 (which in practice 
generates errors ~ 10 -14 ). The radial-transition detection 
routine numerically integrates an object if its osculating or- 
bit as determined on the outer step comes within 5% of the 
transition zone. 

3.1 Single planet 

Here we consider a two-zone integrator with outer timestep 
0.05 yr and inner timestep 0.025 yr. First we set the tran- 
sition zone from 0.9 to 1.1 AU. We place objects of masses 
ranging from 1 M^g to 10000 at semimajor axes between 
0.8 and 1.2 AU, and consider the resulting behaviour of the 
energy (in the one-planet case, a measure of the variation 
in semimajor axis). For objects which enter the transition 
region, the surrogate Hamiltonian under which the planet is 
moving contains terms from both the inner and outer zones. 

Figure [3] shows three points for each of the five mass 
cases (1,10,100,1000, and 10000 M ffi ) at each initial semi- 
major axis: one for the basic DH integrator (eq. I14[) with a 
timestep of 0.025 yr, one for Naoko, and one for the DH in- 
tegrator with a timestep of 0.05 yr. The errors generally de- 
crease with increasing semimajor axis: for a fixed timestep, 
at larger orbital radius there are more steps per orbit. As 
expected, for objects in the transition zone, the new code re- 
ports errors which smoothly interpolate between those of the 
smaller and larger DH runs, and become indistinguishable 
from the standard algorithm outside the transition region. 
In the 10000 Mg^ runs - 3% of the mass of the Sun - there 
are some slight deviations apparent at the edges of the tran- 
sition region. The clean interpolation also breaks down at 
very low mass for a different reason: the integration error is 
dominated by the integration tolerance. 
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Figure 3. Energy error as function of semimajor axis and object 
mass; the DH runs for 0.025 and 0.05 yr timesteps are in green, 
the Naoko runs are in blue. 
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Figure 4. Maximum eccentricity as function of semimajor axis 
and object mass. Most of the 1 M^g runs are not plotted, as their 
maximum e was 0. 



The situation is very similar for the eccentricities, as 
shown in figure [4] Again we see the new code interpolating 
between the 0.05 and 0.025 yr DH runs. This shows that at 
least in a sufficiently smooth case, integrating an isolated 
object on multiple timesteps need not introduce spurious 
energy or angular momentum behaviour. 

Figure [S] shows the relative energy error as a function 
both of the timestep and the inner:outer timestep ratio, 
where the largest outer timestep used was 0.05 yr, for a 
1 Mm object at 1 AU with e = 0.1, with transition zone 
from 0.95-1.05 AU. As expected, the integrator behaves as 
a second-order algorithm. Increasing the ratio to take more 
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Figure 5. Maximum relative energy error as function of timestep 
and innenouter timestep ratio for a 1 object at 1 AU with 
e = 0.1, with transition zone from 0.95-1.05 AU. 



Figure 6. Maximum relative energy error as a function of tran- 
sition width and eccentricity; the eccentricity increases between 
curves by a factor of 1.48. 



inner steps per outer step improves the energy error very 
little beyond 4:1, as the maximum error is controlled by the 
outer step which is not changing. 

Although this condition is necessary, it is clearly insuf- 
ficient. More realistic tests involve objects which repeatedly 
cross from the inner zone to the outer zone through the 
transition zone during the same orbit. 

Without loss of generality we placed a 1 object in 
the centre of a transition zone at 1 AU, and varied both the 
eccentricity e and the relative transition zone half-width h, 
where h is defined so that the zone extends from a(l — h) to 
o(l + /i), in analogy with the eccentricity. The outer timestep 
was chosen to be l/20th of the orbital period, the eccentrici- 
ties were varied from 0.001 to 0.37 (the latter chosen because 
that corresponds to a perihelion for which 0.025 yr is l/20th 
of the orbital period at that distance), and the half- width 
was varied from 0.01 to 0.50 AU. The resulting maximum 
relative energy error is shown in figure [6] as a function of the 
e/h ratio; the contours correspond to constant e. 

At small e/h- that is, when the transition zone is larger 
than or comparable to the radial excursion per orbit - the 
energy error is well-behaved for a fixed e, and is only weakly 
dependent on h. At e/h > 10, there are many integrations 
which do not converge. Nonconvergence is easily recognized 
from the evolution of the semimajor axis or the eccentricity: 
the planet begins to migrate away from its original location 
until it finds a semistable configuration. Examples are given 
in figure [7] for e ~ 0.25. Not only does the transition need 
to be sufficiently smooth along an orbit, but for a run with 
a larger tolerance or a larger timestep, a wider transition 
zone may be required so that the numerical integrator can 
detect the transition (unless the code accepts hints about 
the location of difficulties in the integrand). 




h = 0.031 AU 



h = 0.053 AU 




10 20 30 40 10 20 30 40 

time [kyr] 

Figure 7. Relative change in semimajor axis and eccentricity for 
various transition half-widths h at e ~ 0.25; the three half-widths 
below did not converge. 



3.1.1 Varying the switch 

We also explored the effects of alternate choices of transi- 
tion function. We place a 1 Mgg object at 1 AU with e = 0.2, 
and a transition zone extending from 0.95 to 1.05 AU with, 
with outer timestep 0.05 yr and inner timestep 0.025 yr. 
We then vary the switch over the six functions described in 
£ 12.61 the step function; the switch f(x) = x; the DLL98 cu- 
bic; the DLL98 septic; the Chambers switch; and the Rauch 
& Holman switch. Figure [8] shows the resulting relative en- 
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Figure 8. Maximum relative energy error for single planet at 1 
AU, e=0.2, with transition zone 0.95-1.05 AU, for various switch 
functions. 



ergy errors: the behaviours for the convergent integrations 
continue for millions of orbits. 

The step function fails immediately: the object migrates 
outward in semimajor axis to ~ 1.00456 AU over the first 
2700 years and then stops, corresponding to a maximum 
relative energy error of ~ 0.0046. The remaining integrations 
all succeeded, with /dll3 showing the smallest energy error 
and /rh the largest. We have found this is typical. That 
the crude /x performs better than the smooth /rh can be 
understood by returning to figure [1] /rh is the smoothest 
function, but for a fixed transition width it is also the one 
with the narrowest effective transition (extending roughly 
from x=0.25 to x=0.75). 

The execution times for the various switches (excluding 
the step function) were generally comparable - mostly agree- 
ing within the scatter - except that the Rauch & Holman 
switch was consistently the slowest, and fx the fastest. The 
slower speed of the RH switch is likely chiefly due to the 
evaluation of the trigonometric function. (Note that some 
versions of the common GCC C compiler will occasionally 
refuse to inline functions for obscure reasons, leading to 
strange profiling results.) 

Thus we chose /dll3 as the best compromise between 
efficiency of evaluation and energy conservation. 

3.2 Multiple planets 

The above results are unsurprising, insofar as the use of 
smooth transitions to move portions of the Hamiltonian 
from one part of the step to another is familiar. Of more 
concern is whether the reduction in the number of force eval- 
uations will lead to unacceptable resolution of the angular 
momentum exchange between planets in different zones. We 
will concentrate on the interactions which are likely to be 



the most sensitive to changes in sampling frequency, namely 
resonant interactions. Recall that we are willing to accept a 
cruder approximation to the dynamics than we would ordi- 
narily do, as long as there are no clear integration failures. 
Since one of our interests is studying migration scenarios, in 
which resonances play a significant role, it is important to 
verify that such behaviour is not lost when using the multi- 
zone methods. 

We will verify that resonant behaviour can persist when 
the objects are in different zones, when one object is in a 
transition zone, and when the objects migrate across a zone 
in resonance (coorbital or otherwise). In each simulation we 
have used the two-zone 0.05/0.025 yr 0.9-1.1 AU transition. 
To make comparisons with the behaviour of the DH inte- 
grator clearer we have disabled the encounter treatment for 
the following sections (until t|3.2.7|l . 

3.2.1 2:1 mean motion resonance 

We place two 2 Mgg objects on cold e ~ 0.002, i ~ orbits 
in the 2:1 mean motion resonance in three configurations: 
(1) ai = 0.63 AU, <J2 = 1.00 AU, where the outermost ob- 
ject is in the middle of the transition region; (2) ai = 0.76 
AU, (22 = 1.20 AU, where neither object is in the tran- 
sition zone but the objects are in different zones; and (3) 
ai = 1.00 AU, a2 — 1.59 AU, where the innermost object 
is in the transition zone. (All semimajor axes are approxi- 
mate; the initial configurations were found by introducing 
a slowly-decreasing dissipation into the DH algorithm.) The 
evolution of the resonance angles is plotted in figure [5] In 
all cases the resonance is preserved, and in the cases where 
one object was in a transition zone the agreement is excel- 
lent between the Naoko results and the results of a DH run 
with timestep 0.025 yr. In the second case - in which neither 
object is in a transition zone - the libration is considerably 
larger in the Naoko run. This is not surprising: in the ab- 
sence of dissipation the resonance is quite sensitive to the 
initial conditions, and case 2 has the greatest sudden change 
from DH to Naoko, as none of the force is being evaluated 
on the higher (inner) frequency. Even in this case, the in- 
troduced libration is comparable to the differences in width 
between different DH runs in which only the initial angles 
had been changed. For our purposes it is sufficient that the 
objects remain resonant, which they do. 

3.2.2 Co-orbitals 

We place two objects - one 10 and one 1 - at 1.3 
AU on (osculating) circular and coplanar objects with mean 
anomalies differing by 1.5 radians, and let them migrate in- 
wards in the standard Hayashi minim um mass solar nebula 
(MMSN) disc model (jHavashil 11981m with surface density 
~ 1700 g/cm 2 (r/AU) -1 ' 5 .) The migr ation was driv e n by a 
prescription for type I gas migration iTanaka et al.l |2002), 
in which the motion is produced by the asymmetry between 
the torque s generated by the interior and exterior w akes of 
the body l|Goldreich fc Tremainel Il980l : IWard||l986h . (The 
difference in mass between the two test objects is to ensure 
that the lower-mass planet is being carried along by the 
resonance and not merely migrating in tandem.) Figure [101 
shows the evolution of the semimajor axis and the resonant 
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Figure 9. Evolution of resonance angles for pairs of objects in 
2:1 mean motion resonance under both DH (green) and Naoko 
(blue) integrators. The three cases correspond to a\,a^ = (0.63 
AU, 1.00 AU), (0.76 AU, 1.20 AU), and (1.00 AU, 1.59 AU). 



Figure 11. Scmimajor axis and resonance angle for pair of ob- 
jects in 3:2 resonance crossing transition zone (the shaded area) 
for DH (green) and Naoko (blue) integrators. 




time [yr] 



Figure 10. Semimajor axis and resonance angle for coorbital 
migration pair crossing transition zone (the shaded area) for DH 
(green) and Naoko (blue) integrators. 



angle for both the DH and Naoko integrators. Even when the 
pair crosses the transition zone, the trajectories are barely 
distinguishable. It may seem counterintuitive that the algo- 
rithm performs near-perfectly on coorbital resonances which 
depend sensitively on the local potential and considerably 
worse on distant resonances which might be expected to be 
more forgiving of sampling errors. However, the method's 
errors compared with the the standard DH algorithm all in- 
volve relative differences, and in the case of co-orbital objects 
on near-circular objects, these are minimal. 



3.2.3 Transition crossing 

A more rigorous test comes from seeing whether resonance 
can be maintained when the objects are under different - 
and varying - terms in the Hamiltonian. We place a 10 
planet at 1.40 AU and a 20 M ffi planet at 1.85 AU, just 
outside the 3:2 resonance at ~ 1.83AU, and turn on type 
I migration in the MMSN for them both. The results are 
shown in fig. 1111 The 20 planet migrates faster than the 
10 planet and catches up to the resonance within the 
first ~ 1000 yrs. The evolution under DH of 0.025 yr and 
Naoko with 0.025/0.050 yr is almost identical until the in- 
ner object enters the transition zone; that Naoko is using 
a timestep twice as large makes little difference. Even after 
the objects enter the transition region, the planets remain 
in the 3:2, and the libration width is comparable between 
the two runs. 



3.2.4 Transition crossing: varying timestep ratios 

It is of interest to see how the maintenance of the resonance 
varies as the inner:outer timestep ratio is varied. We place 
terms in the Hamiltonian. We place a 1 M(g planet at 0.77 
AU and a 10 planet at 1.23 AU, just outside the 2:1 
resonance, and again turn on MMSN type I migration. We 
fix the outer timestep at 0.05 yr, and change the number of 
inner steps per outer step, while keeping the transition zone 
from 0.90 AU to 1.10 AU in each run. Figure [12] shows the 
resulting evolution of the resonance angle. The use of mul- 
tiple timescales for the forces introduces spurious libration 
spikes (with corresponding behaviour in the eccentricities), 
with the libration width during the spikes generally increas- 
ing with the timestep ratio. However, in all cases, from 2:1 
through 96:1, the integration error failed to break the res- 
onance during the transition, and the pair returned to the 
correct trajectory after the outer object passed through the 
transition region. 
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Figure 12. Evolution of resonance angle for 2:1 pair crossing a 
transition region as the inner:outer ratio is varied; the thin line is 
the reference single-timestep behaviour 



3.2.5 Transition crossing: effect of force frequency 



As mentioned in £12.51 although we have chosen to evaluate 
the forces on the outermost timestep, one can vary this fre- 
quency. We repeat the test of the previous section but with 
the forces always evaluated on the inner timestep (and ac- 
cordingly with no need for any transition function applied 
to the kick operators), in which the only active transitions 
involve the timestep on which the object is being drifted. 
The resulting evolution of the resonance angle is shown in 
figure [TBI We see that using the fixed high-frequency forces 
results in a considerable qualitative improvement, although 
the multiple timescales used for the Kepler drift still gen- 
erate libration spikes. This serves to bound the likely im- 
provement we can imagine by increasing the force frequency, 
which is very costly. Recall that the formation scenarios we 
seek to study using the new method involve large particle 
numbers where the execution time is dominated by the force 
calculation - a quite different regime from the solar system 
modelling of ISaha fe Tremaind l| 1994) . 



3.2.6 Convergence 

It is important to recognize that by changing the time sam- 
pling as we have done, we are unavoidably changing the 
dynamics of the system. Since the outer timestep remains 
fixed as we increase the inner timestep, and with it some 
fraction of the force evaluation, we are modifying the ef- 
fective interaction between the inner and outer objects. For 
example, we consider a pair of 1 M^g objects near a nomi- 
nal 2:1 resonance at ~0.79 AU and 1.25 AU and vary the 
timestep and innenouter ratio. The original outer timestep 
was 0.05 yr. Figure [14] shows the resulting change in eccen- 
tricity evolution for the inner object as both the inner:outer 
timestep ratio and the outer timestep are varied (the outer 
object shows similar behaviour.) As the ratio is increased 
in the top section of the figure - i.e. as the inner timestep 



Figure 13. Evolution of resonance angle for 2:1 pair crossing for 
different transition region as the inner:outer ratio is varied; the 
heavier line corresponds to the high frequency force evaluation, 
and the lighter line to the standard low frequency evaluation of 
figure [12] 



is decreased - the eccentricity evolution does converge, but 
it does not converge to the correct trajectory, although the 
relative differences in this case are not large. For a fixed 
timestep ratio, of course, as the outer timestep is decreased, 
the evolution converges as it should. 



3.2.7 Close encounters 

As mentioned previously, when two objects are in the same 
zone the close encounters are treated as in the Chambers al- 
gorithm, and we disallow encounters between objects which 
are not in neighbouring zones. Since it is known that ac- 
curate resolution of close encounters is quite sensitive to 
the effective Hamiltonian under which it is being integrated 
(DLL98, §4), it is reasonable to expect poor behaviour for 
objects undergoing close encounters in a transition region 
where it is undergoing drifts of different duration. A de- 
crease in accuracy may be tolerable at these few locations 
in a simulation, but the objects cannot suffer dramatic in- 
stabilities. 

We use the 0.9-1.1 AU, dt=0.05 yr/0.025 yr transition. 
To generate frequently-encountering systems, we placed a 5 
M0 planet at 1.00 AU, in the middle of the transition re- 
gion, on a circular orbit and a protoplanet of 0.1 Mgg at 1.5 
AU with e = 0.333 on a coplanar orbit, set all angles to 
and varied only the mean anomaly from to 27r. After re- 
moving runs in which the smaller body did not undergo an 
encounter, would have suffered a potential merger, or (in one 
case) became sufficiently eccentric that it escaped the inte- 
gration region, 16 runs remained. The relative change in the 
energy and the Tisserand p arameter T = l/2a+ \J a(l — e 2 ) 
jMurrav fc Dermott][l999T ^3.4) is plotted in figure [151 The 
RMS of the maximum changes in E is ~ 0.0004, and that 
of T ~0.01; the horizontal segments correspond to periods 
without encounters, where the conservation properties re- 
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Figure 15. Energy conservation and variation of Tisserand pa- 
rameter during close encounters described in i|3.2.7l 



Figure 14. Eccentricity evolution for inner object in near- 
resonant pair described in £13.2,61 for various timesteps and in- 
neriouter timestep ratios; the system does not converge to the 
true path as the ratio is increased but does as the timestep is 
decreased. 



vert to the standard non-encountering behaviour. For com- 
parison, the equivalent RMS maximum changes for our refer- 
ence SyMBA implementation with timestep 0.025 yr - with 
mergers disallowed, to make the integration even more diffi- 
cult - were ~ 3 ■ 10" 6 and ~ 0.003 for E and T, respectively. 
It is clear that the repeated close encounters in the transition 
zone cause a considerable decrease in accuracy, especially in 
the energy, but the error growth should be tolerable. 

The Bulirsch-Stoer numerical integrator we use often 
fails, and reports its failure, at high precision for very close 
encounters. This is typically not a problem in practice, as 
the separations at which the integrator fails are usually much 
smaller than the distances at which we would merge the two 
bodies, but is an issue for small-radius dynamical problems. 
The same failure mode exists when no transition zones are 
involved and the symplectic algorithm is one resembling that 
of Chambers. 



3.3 Full example: short-period Neptunes 

As a case study, we consider the problem of forming short- 
period Neptunes, which are particularly mysterious. It is 
clear they did not form at the small orbital radii at which 
they are currently found, and so it is likely they have 
migrated in from further out. However, their masses - 
1OM0 and higher - are unexpected. They are large enough 
that the timescale for type I migration to drag them into the 
Sun is considerably shorter than the formation timescale, 
but too small for some plausible mechanisms for suppress- 
ing type I migration to be effective (such as the opening of 
a gap in the disc at the transition to type II migration.) 
There are two main scenarios present in the literature 



for fo rming these objects. The first, due to lAlibert et al.l 
(2006), concentrates on sophisticated gas physics and follows 
the evolution of effectively isolated cores through the disc 
as they grow via planetesimal accretion and migrate, and 
finally have their atmospheric mass reduced after arriving in 
the short-period region by eva poration. The second, due to 
iTerquem fc Papaloizoul |2007), takes a more traditional N- 
body approach, and follows the evolution of a small number 
of large embryos in the inner region. Neither of these are easy 
to reconcile wit h the standard oliga rchic growth model for 
core formation (K okubo fc Idall 19981 ): the first requires core 
formation to be suppressed everywhere in the disc except at 
a few special locations, and the second uses initial conditions 
for the mass distribution which it is difficult to recover self- 
consistently from earlier stages of oligarchic evolution. 

It is of considerable interest whether the standard model 
can generate anything resembling the observed planets - if 
it cannot, then we have strong evidence that we are missing 
important physics - and the new method can address this 
question. We will merely sketch the application here; forth- 
coming work will describe our results in more det ail. We 
apply the two-stage approach of McNeil et all l|2005f ). which 
was based on that oflTh ommes et all (|2003l )T using semiana- 
lytic approximations of oligarchic evolution to treat the early 
evolution of the disc and then generating an N-body real- 
ization when the number of N-body particles needed drops 
to practical levels. We consider various disc models (a ll re- 
sembling the minimum mass model of lHavashil Il98lh and 
different migration efficiencies. 

An example of a low-resolution run is shown in figure 
[16] which employed 21 embryos and 212 planetesimals. The 
semianalytic model was run to 1 Myr for a disc with sur- 
face density proportional to r -0 ' 5 of ~3 times the standard 
Hayashi mass without a snow line, with nominal type I mi- 
gration and aerodynamic drag (for a 1 km planetesimal size), 
where the e-folding time for the decay of the gas disc was 
0.5 Myr. A section of the system ranging from 1.0 AU to 
5.1 AU was then evolved under the new N-body code until 
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the gas was mostly absent. As a consequence of the semian- 
alytic model predicting that the innermost part of the disc 
had gone to oligarchic completion, the embryos inside of 2.0 
AU migrate inwards smoothly in tandem: there is nothing 
to perturb them. This results in 10.3 of embryo mate- 
rial inside of 1 AU when the gas vanishes in an apparently 
stable configuration. Further out, the system is still undergo- 
ing chaotic evolution. In a more realistic situation, in which 
there were still planetesimals around in the inner regions to 
perturb the embryos, they might merge, producing larger 
planets; or, if they merge too early, they could migrate into 
the Sun. 

Figure[T7]presents examples of the kinds of interior con- 
figurations our toy simulations produce; all are plotted at 
t=3.25 Myr. Typical results include failures in which no or 
almost no mass is left in the interior region (S2, S5, S7, 
S9); cases where the total mass is interesting from the per- 
spective of hot Neptune formation but the embryo mass is 
spread out over a large number of embryos instead of be- 
ing concentrated in one or two objects (S4, S8); and cases 
where the total mass is appealing and there are only two 
embryos of multiple Earth masses (SI, S3, S6) but where 
too much of the mass is locked up in planetesimals, possi- 
bly as a consequence of the low resolution. (Planetesimals 
do not self-interact in these simulations.) Even these crude 
simulations are useful to estimate whether the conditions 
are appropriate for larger integrations. We will present re- 
sults from higher resolution simulations in a forthcoming 
publication. 



4 DISCUSSION 

The observed behaviour of the method (including on tests 
not presented above) is compatible with the principles that 
(1) parts of the Hamiltonian can be moved between op- 
erators with different timesteps if the transition is suffi- 
ciently smooth, (2) reversibility is more important than ac- 
curacy for secular energy conservation, and (3) decreased 
time resolution is not catastrophic for dynamics on much 
longer timescales. We emphasize that we do not claim, and 
have not shown, that the integrators are in any sense op- 
timal. (Some evidence suggests that a drift-kick-drift split- 
ting can be preferable to the kick-drift-ki ck splitting used 
here, for example; see I Wisdom et al. I ll996l .) To use a loose 
analogy from the history of close encounter treatment, this 
appro ach probably lies somewhe re between the hybrid meth- 
ods of lLevison fc Duncan! (| 1994 ) and the mature methods of 
DLL98. As always, one must bear the approximations be- 
ing used in mind - small effects can accumulate to cause 
dramatic changes on secular timescales in ways not always 
easy to predict. The approximations used here have the 
virtue of being different from those used previously. Given 
the decrease in force accuracy involved, it may be best to 
think of the use of these methods as one would a reversible, 
Kepler-adapted treecode: useful for studying certain for- 
mation problems, and generally informative in a statistical 
sense, but inappropriate for cases where detailed dynamics 
are important. Potential users are advised to think carefully 
about whether the approximations are suited for their prob- 
lems, as there are far more ways for an integration to fail in 
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Figure 16. Embryo evolution with time for a simulation of a disc 
3 times the minimum mass with £ <x r~ °' 5 , with disc dissipation 
c-folding time of 0.5 Myr. For each embryo lines corresponding 
to semimajor axis, perihelion distance, and aphelion distance are 
drawn. 



the multi-zone case than there are in standard SyMBA-style 
integrators. 

There are several potential applications and directions 
for generalization we have not considered here. As dis- 
cussed in section |2"31 iLevison fc Duncan! (|2000h modified the 
Chambers approach to handle objects which occasionally re- 
quire very small timesteps, at the price of integrating the 
entire system numerically when such an approach occurs. 
The new method succeeds in decoupling such objects, and 
so may be practical in some cases for which their method is 
inapplicable (such as when ejections of low-mass objects are 
common). It may also be useful in studying planet forma- 
tion in close binaries, for which some sympl ectic methods for 
deali ng with these systems already exist [[Chambers et al] 
120021). 

ISaha fc Tremainel (|l994l ) develop 'symplectic interpola- 
tion' methods to improve the accuracy of their integrator, by 
using a symplectic prediction of the effects of Hkc P to bet- 
ter synchronize the force calculations. This would probably 
work here as well, in a formal sense, but it is unlikely to pro- 
vide dramatic improvements, and the errors in energy and 
angular momentum conservation introduced by the method 
should be tolerable for most formation problems as they 
stand. 
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Figure 17. Mass versus semimajor axis snapshots at T=3.25 
Myr for various runs; the larger blue circles represent embryos, 
and the smaller red circles represent planetesimals. The total mass 
in embryos (Mt) and planetesimals (raj) in the region inside 1 
AU is given in the upper left. 



5 CONCLUSIONS 

We have presented a new integration method which pre- 
serves most of the speed and conservation properties of stan- 
dard close-encountering symplectic integrators for planetary 
dynamics, but introduces radial zones between which the 
drift and force evaluation timesteps can vary. This allows for 
a new tradeoff between the accuracy of the integration and 
the speed, one which is inappropriate for precise dynamical 
studies but may be of use for investigating other scenarios. 
We expect it to be useful for approximate N-body modelling 
of planetary formation problems where (1) there is a wide 
range of orbital velocities, (2) the majority of the material 
needing particle representation is in the regions with longer 
orbital periods, and (3) the system remains sufficiently cold 
that the transition regions between radial zones can be kept 
reasonably small. Migration scenarios of the formation of hot 
exoplanets satisfy all these constraints, and we are currently 
exploring such applications. 
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